EPJ manuscript No. 

(will be inserted by the editor) 



Matrix formalism and singular value decomposition for the 
location of gamma interactions in segmented HPGe detectors 

p. Desesquelles\ T.M.H. Ha\ K. Hauscliild\ A. Korichi^ F. Le Blanc^ A. Lopez-Martens\ A. Olariu^ and CM. 
Petrache^ 

^ CSNSM CNRS/IN2P3 and Universite Paris Sud 11, 15 rue G. Clemenceau, 91405 Orsay, France 
2 IPNO CNRS/IN2P3 and Universite Paris Sud 11, 15 rue G. Clemenceau, 91406 Orsay, France 

On behalf of the AGATA Collaboration 

Received: date / Revised version: date 

Abstract. Modern coaxial and planar HPGe detectors allow a precise determination of the energies and 
trajectories of the impinging gamma-rays. This entails the location of the gamma interactions inside the 
crystal from the shape of the delivered signals. This paper reviews the state of the art of the analysis of the 
HPGe response function and proposes methods that lead to optimum signal decomposition. The generic 
matrix method allows fast location of the interactions even when the induced signals strongly overlap. 
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gamma-ray spectroscopy - 07.50.Qx Signal processing electronics 



1 Introduction 

Modern segmented germanium detectors have become an 
indispensable tool for both the measurement of gamma- 
rays energies [T.T.'S and the 3D positioning of gamma-ray 
sources. Germanium detectors have evolved a lot since the 
sixties, notably by the increase of the crystal size and the 
possibility to cut the crystal into electrical segments. As 
for many physical detection devices [4,5 ,6,7 , the present 
improvements are mostly due to the development of on- 
line pulse shape analysis techniques. 

The segmentation of the germanium crystal is induced 
by electric contacts distributed on its surface. The volume 
of each resulting electrical segment is of the order of a 
several cubic centimeters. 

Two main classes of such segmented germanium detec- 
tors exist. The first class, coaxial detectors, allow a precise 
determination of the energy (up to several MeV, for crystal 
length of the order of 10 cm) of the gamma-rays emitted 
by a nucleus and, in the case of a moving nucleus, their 
angle of emission so that measured energies may be prop- 
erly Doppler corrected. In this frame, the detector array 
plays the role of a theodolite and a calorimeter. It must be 
able to detect in coincidence a large number of gamma- 
rays and to measure the characteristics of each photon 
individually. 

The arrays of planar or coplanar-segmented detectors 
play also the role of a positioning system of gamma sources. 
The ability of segmented germanium detectors to locate 
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gamma-ray sources makes them valuable instruments for 
medical imaging, radioactive source search in the frame 
of national security or environmental monitoring. Planar 
array are typically made of layers of parallelepipedic seg- 
ments. From the first two interactions of the incident gamma- 
ray, the Compton formula permits to determine a cone of 
possible directions to the source. The intersection of many 
such cones defines the position of the source [8l[9]. 

The location of the gamma-ray interactions (hits) in- 
side the crystal is performed using pulse shape analysis. 
Indeed, Compton scattering, pair creation or photoelectric 
absorption generate a number of electron/hole pairs. The 
created charges migrate in the electric field to the contacts 
at the surface of the crystal. This motion induces a vary- 
ing image charge on the electrodes belonging to the seg- 
ment where the hit occurred and on the electrodes of the 
neighboring segments. The amplitudes of the signals are 
proportional to the deposited energy and the amplitude 
of the signals induced on the cathodes of the neighboring 
segments increases with the proximity of the hit. Hence, 
the location of the hit and the energy deposit can be de- 
duced from the shape of the signals. When more than one 
hit occur simultaneously in the same segment or in neigh- 
boring segments, the total signal is the sum of the individ- 
ual signals. Thus, the location of the interactions can be 
determined using a signal decomposition algorithm. The 
mathematical bases of such algorithms will be discussed in 
the following. Signal decomposition appears to be a more 
difficult problem in the case of coaxial detectors then in 
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the case of planar detectors due to the variety of segment 
shapes and the complexity of the electrical field map. 

When the locations and the energy deposits have been 
determined for all hits, a so-called tracking code p^O l fTT ] 
[T2] rebuilds the paths of the gamma-rays from hit to hit. 
Finally, the sum of the energy-deposits and the location of 
the first hits gives the characteristics of the gamma-rays. 

This paper is organized in the following way. In the 
next section, we introduce the matrix formalism which 
allows to describe the mathematical link between signal 
shapes and locations of the hits. In the next section, the 
properties of the response function of the germanium de- 
tectors is deduced from the analysis of the transformation 
matrix. It will be shown that signal decomposition is an 
ill-posed problem. Moreover, signal decomposition entails 
the solving of a very large set of linear equations. The 
Singular Value Decomposition (SVD) method permits to 
solve both problems and to speed up the decomposition 
codes. In the forth section, we describe a complete proto- 
col for on-line signal decomposition. 
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Fig. 1. Example of meta-signal. The fourth signal corresponds 
to the hit segment and the eight other signals correspond to 
its neighbors. 



2.2 Linear system of equations 



2 Mathematical formalization 



Using the linear and additive properties of meta-signals, 
one obtains: 



2.1 Introduction 

Our goal is to determine the locations (xi^yi^Zi) and the 
energy deposits of the M interactions of the gamma- 
ray(s) inside the germanium crystal, knowing the sum 
of the M individual signals. This signal decomposition is 
made possible thanks to two properties: the amplitudes of 
the individual signals are proportional to the correspond- 
ing energy deposits, and the signals are additive, that is 
the signal actually delivered by the segments is the sum of 
the signals induced by each gamma interaction. A given 
hit induces signals in the segment where it occurred and 
in the neighboring segments. Thus it is useful to introduce 
the notion of meta-signal as simply the concatenation of 
the hit segment signal and its neighbor segment signals. 
An example is shown in fig. [l] where the signal of the hit 
segment runs from samples 157 to 208 and is concatenated 
with the signals of its eight neighbors. The number of sam- 
ples (52) is chosen so that it includes the rise time of the 
hit segment signals, which maximum value is 37 samples, 
and at least 15 samples at the minimum of the signal as 
it gives the deposit energy. The number of neighboring 
segments is five for the first an the last layers of coaxial 
detectors. 

The meta-signals will be denoted s{t). The meta-signal 
can also be seen as a vector s whose components are simply 
the amplitude of the signal in the successive bins Q This 
vector represents the whole information delivered by the 
detector. In the following, the meta-signal corresponding 
to a unit energy deposit at point (x, z) will be noted 
m{x,y,z,t). 



^ In the following, vectors will be denoted by small bold let- 
ters and matrices by capital bold letters. 



s{t) = JJJ y, z) m{x, y, z, t) dx dy dz , (1) 

where e{x^y,z) is the energy deposited at point {x^y^z) 
and the sums run on the volume of the detector. If the 
meta-signals m{x^y^ z^t) are known, then the inversion 
of eq. gives the energy deposits and their locations. 
The task of signal decomposition is to solve this so-called 
inverse problem. In actual applications, the meta-signals 
are calculated using a simulation code [13] or a crystal 
scanning system [T4l[T5]. Thus, they are known only for 
discrete points {xj^yj^Zj) on the nodes of a given grid. 
The step of the grid is typically of a few millimeters. The 
previous equation becomes: 

s{t) = ^ej m{j,t), (2) 

j 

where ej is the energy deposited in the voxel surrounding 
the grid point j as shown in fig. |2j 

Noting nij the basis meta-signal corresponding to a 
unit energy deposit at point j, this equation can be rewrit- 
ten as: 

j 

and finally in a matrix form as: 

Me = s, (4) 

where M is the transformation matrix. The j*^ column of 
M is m^- (fig.[3|. 

Signal decomposition consists in solving this matrix 
linear system, that is to find the components of the e 
vector. An example of energy vector is shown in fig. [2] 
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Fig. 2. Example of calculation of the hit estimated location 
in a 2D rectangular segment. The open circles stand for the 
grid points. The dotted square around grid point 6 represents 
its voxel. The star indicates the actual location of the hit. The 
surfaces of the dots are proportional to the corresponding com- 
ponent in the e vector. The energy weighted bary center (esti- 
mated location of the hit) is indicated by a diamond. The error 
of the method (resolution) is the distance between the diamond 
and the star. 



M 



Fig. 3. Graphical sketch illustrating eq. (|4|. The bold box 
stands for the M matrix and the vertical segments for the 
e and s vectors. The logic of such diagrams is explained in 
Appendix A. 



In principal, the general solution of eq. Q is obtained 
by inverting matrix M [18 . Yet, the mere inversion method 
would require a very long computation time and would 
give an unphysical solution. In the following, we will see 
how these difficulties can be overcome. 



2.3 Properties of the solution vector 

The components of the solution vector represent the en- 
ergy deposit in each voxel of the grid. Of course, most of 
the Cj components are equal to zero, and, as they repre- 
sent energies, the non-null components are positive q In 
the hypothetical case when the hit occurs exactly on a 
grid point, only the corresponding solution component is 
different from zero and, as the basis signals correspond to 
unit energies, this component is equal to the energy de- 
posit. In the real case, the hit location does not belong to 
the grid. The strategy of most algorithms, such as the grid 
search, is to find the signal from the grid that best matches 
the measured signal. The corresponding grid point is con- 
sidered to be the closest to the actual gamma hit. In this 
case, the precision on the hit location is directly connected 
to the grid step. Using the matrix formalism, a more ac- 
curate location is possible. Indeed, after solving of eq. Q, 
more than one component of the solution vector can be 
different from zero, even for a single hit. The estimator of 
the hit location is then calculated as the energy weighted 
barycenter (fig. [2| of the components: 



(5) 



For this purpose, many algorithms have been devel- 
oped, (grid search p^ fTT], matrix method p ^ [T9 l l2Q], wavelet 
decomposition [21], ...) some of them using artificial intel- 
ligence methods (genetic algorithms [22^, neural network 
[23] , ...). In fact, artificial intelligence methods appeared 
to be efficient but slow and thus applicable only in the 
case when the interaction location is performed off-line. 
For real-time applications, the most used methods are im- 
proved forms of grid search. Unfortunately, these methods 
are not well adapted to the situations when more than one 
hit occurs in the same segment or in neighboring segments. 



where the Xj are the three dimensional locations of the 
grid points. The denominator is equal to the deposited 
energy. 

Of course, due to the uncertainties on the measured 
signal (electric noise, time alignment [23][24] |25] , cross-talk 
[26 , etc.) and on the transformation matrix (grid dis- 
cretization), it is not possible to find the exact position 
of hits. Both the resolution (i.e. the r.m.s. distance be- 
tween the estimated and the actual hits) and the spread 
of the cloud of non-zero components of the solution vector, 
shown as dots in fig. [2] depend on the alterations of the de- 
tected signal. As an example, fig. [4] shows the evolution of 
the resolution and of the cloud spread as a function of the 
signal-to-noise ratio and of the time jitter of the signals. 
The time jitter applies to the met a signal as a whole. The 
time jitter between segments is much smaller and does not 
influence the resolution. The data are simulated using the 
MGS code. The deterioration of the resolution appears to 
be mainly due to the noise whereas the spread of the cloud 
is also sensitive to the amplitude of the time jitter. Large 
cloud spreads do not allow the separation of hits, which 
lie close together, as will be seen in section 4. 

In the following, we focus on the mathematical prop- 
erties of eq. Q and on its optimum solving method. 



^ The total energy deposited in a segment is directly deduced 
from the amplitude of the charge signal. 
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Fig. 4. Resolution (top panel) and spread of the cloud of non- 
zero components of the e vector (bottom panel) as a function 
of the signal-to-noise ratio and of the standard deviation of the 
signal time jitter. These results are obtained by the NNLS solv- 
ing of eq. Q using MGS simulations of AG ATA signals. As- 
suming a noise level of 3 keV sigma, the signal-to-noise bins cor- 
respond, from left to right, to about 10 keV, 33 keV, 100 keV, 
333 keV and 1 MeV gamma hits. 



2.4 Least square solving 

The number of unknowns in eq. Q (columns in the M 
matrix) is the number of grid points and the number of 
equations (lines in the M matrix) is the number of samples 
in the meta-signal. This consideration gives a first upper 
limit to the number of grid points that can be accommo- 
dated so that the linear system is not underdetermined. 
Hence, the M matrix for a segment can be made vertical 
rectangular. Due to the signal alterations, eq. (4| has no 
exact solution, thus, we are searching for the solution egoi 
that minimizes the residue ||Mesoi — s||. 

If there were no constraints on the components of the 
solution vector then the least square solution would be 
given by: 



Gsoi = CMM)--' *M s 



(6) 



where *M is the transpose of matrix M. In fact, this al- 
gebraic solving is not valid since, due to the noise and the 
uncertainties that affect the system, the resulting solution 
would most probably have negative components. Differ- 
ent algorithms are meant to find non-negative solutions 
to this kind of linear system (such as backtracing [27 ) 
but most of them are too slow (the signal decomposition 
must be realized on-line, that is in a few milliseconds). To 



our knowledge, the best compromise algorithms are NNLS 
(nonnegative least squares [28 ) and NNLC (nonnegative 
least chi-square [29 ), as they maximize the number of null 
components and minimize the size of the matrices to be 
inverted and thus the computing time. 



2.5 Uniqueness of the solution 

We will now discuss the characteristics of the matrix M 
that entail the uniqueness and the stability of the solu- 
tion. The properties of M are induced both by the re- 
sponse function of the segments and by the choice of the 
number and the locations of the grid points. When the 
response function is not bijective then, whatever the algo- 
rithm and grid choices, the hits cannot be unambiguously 
located. This situation arises in some particular cases in 
coaxial and large segment planar detectors: the signal re- 
sulting from the addition of two hits may be very simi- 
lar to the signal resulting from a single interaction at the 
barycenter of the two hits (see fig. |5|. This explains why 
the determination of the number of hits within one seg- 
ment is often difficult and some times impossible. If one, 
however, is searching for single hit event^ then, what- 
ever the found solution, the barycenter of the components 
is approximately the same (fig. Isl) . 

The fact that the response function does not have a one 
to one relation between resulting pulse shape and grid lo- 
cations, is equivalent to the fact that, whatever the chosen 
grid size, the matrix M is ill-conditioned. This does not 
mean that more than one solution is the usual situation. 
Indeed, on the one hand, the number of solutions is re- 
duced by taking into account the physical constraints: the 
components have to be non-negative and most of them are 
equal to zero. On the other hand, indiscernability occurs 
only in the low sensitivity zones of the detector. 



3 Solution using singular value decomposition 
3.1 SVD principles 

We have seen that the M matrix has two main defects: on 
the one hand, it is large, thus the calculation of its inverse 
is irretrievably long and, on the other hand, small fluctua- 
tions on the detected signal induce large uncertainties on 
the hit locations. Both of these problems can be corrected 
using singular value decomposition. In the following, we 
will see how this technique can be adapted to our purpose. 

Any square or rectangular matrix can be decomposed, 
as shown in fig. [6) into the product of three matrices: 



M = UW'V, 



(7) 



such that W is a diagonal positive matrix whose diagonal 
components are arranged in a decreasing order, and U and 
V are column-orthonormal matrices (*UU = *VV = 1). 



^ Events are defined as sets of hits induced by gamma-rays 
entering simultaneously the crystal. 
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Fig. 6. Graphical sketch of the singular value decomposition 
in the case of a vertical M matrix (eq. ([t])). The diagonal line 
in the W matrix symbolizes the set of singular values. 
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Fig. 5. Example of two hits giving the same signal as a single 
hit in the low sensitivity zone of a coaxial detector segment. 
The upper figure shows the location of the grid points (small 
open dots), of the single hit (open blue disk) and the two hits 
(red disks). The surfaces of the disks are proportional to the 
energy deposits. The bary center of the two hits is indicated by 
a diamond. The resulting almost over- imposed meta-signals are 
shown in the lower figure, with the same colors (dotted red line: 
two hits, thin blue line: one hit). These signals were generated 
using the MGS [13^ code (neither the noise nor the electronic 
response function are taken into account here). 



The set of singular values that is the values of the 
diagonal matrix, is unique. Singular value decomposition 
has many applications in Physics ranging from the predic- 
tion of the perturbation growth [30], to Principal Com- 
ponent Analysis [31,32 (inertia analysis in multivariate 
space) or genomic analysis [33 . 

Once the matrix is decomposed, the least square solu- 
tion vector could be very conveniently calculated as: 



esoi=VW-i*Us, (8) 

where is the diagonal matrix of the 1/wi (when a sin- 
gular value is equal to zero, the corresponding component 
in the inverse matrix is also zero). Unfortunately, this so- 
lution does not necessarily respect the physical constraints 
proper to our pulse shape analysis problem. However, as 
we will see in the following, SVD remains very useful for 
signal decomposition. 

The first advantage of SDV is that it allows to trans- 
form the rectangular system of eq. Q into a square one 
and to lower the dimension of the matrix to be inverted in 
the case when the number of samples of the meta-signals 
is greater than the number of grid points in the segment. 
Indeed, from eqs. Q, and ([7|), one obtains: 

*Ve = W-^*Us. (9) 

The matrix R = *U has the effect of reducing the 
size of the signal, as shown in fig. [7| in an optimum way, 
i.e. keeping the whole relevant signal-information. The re- 
duced signal will be denotes s^. 

Finally, a new smaller linear system can be substituted 
to eq. Q: 

*Ve = s^ (10) 

The number of lines in the system is decreased as can 
be seen by comparing fig.[8]to fig.|3] Each column of *V is 
the reduced signal corresponding to a unit energy deposit 
on a given grid point. This matrix plays the same role for 
reduced signals as M for signals. 

In the following, we will see that the size of the system 
can be reduced even more using the so-called SVD trun- 
cation. It will be shown that the maximum reduction is 
connected to the condition number of the matrix. 
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Fig. 7. Reduction of the size of the signal with matrix R. 



Fig. 8. Graphical sketch of eq. ([10|. The size of the linear 
system is reduced with respect to fig. |3] 



3.2 Properties of the germanium response function 



where r is the index of the lowest non-zero singular value. 

Two types of ill-conditioned problems exist, which have 
to be addressed with different solving methods. The first 
includes the rank deficient problems which are character- 
ized by a transformation matrix having two, well sepa- 
rated, groups of large and small singular values. In this 
case, the numerical rank is equal to the number of large 
singular values (and r is fixed to this value). The second 
one includes the discrete ill-posed problems, for which the 
set of singular values decreases smoothly. 

In fig. |9j we show the singular values for a 36 segment 
coaxial germanium detector for 2 mm and 5 mm cubic 
grids (here, the M matrices have respectively 41874 and 
2544 columns for 1872 lines). The detector signals are sim- 
ulated using the MGS simulation. The two distributions 
being very close, within a constant factor, we can say that 
the 2 mm grid adds little information with respect to the 
5 mm grid. The condition number of the matrix is of the 
order of 10^^, which is very bad, thus SVD truncation is 
indispensable. The singular values distributions show the 
same drop after one thousand. Thus the value of the cut- 
off r should not be larger than one thousand. In fact, due 
to uncertainties in the nij signals, the cut-off has to be 
even lower, as will be shown in section 3.4. 

A clear gap appears between the first 36 singular values 
and the next one. It means that, from the signal shapes, 
it is very easy to know which segment was hit, and, as 
Wi/wsQ is close to one, this determination is robust with 
respect to the signal noise. This does not come as a sur- 
prise. 



An important characteristic of linear systems is the way 
they transmit the uncertainties, such as the noise, present 
on the right-hand term s, to the solution e. The maximum 
amplification coefficient for the relative uncertainties is 
called the condition number: 

/ l|M-Ms|| \ 

C = max,J^ . (11) 

V ii^ii / 

which is the ratio of the relative fluctuations on the so- 
lution (the fluctuations are measured as the norm of the 
vector) and on the signal. As can be seen in eq. ([8|, the 
amplification increases with the inverse of the singular val- 
ues. In fact, it can be shown that the condition number 
is equal to the ratio between the largest and the lowest 
non-zero singular values. In the case of germanium detec- 
tors, the condition numbers may actually reach very high 
values. Therefore, a prior mathematical treatment of the 
transformation matrix is necessary to ensure the reliabil- 
ity of the solution. Singular Value Decomposition permits 
such a treatment. 

An optimum way to lower the condition number, in- 
troducing a minimum bias on the solution, is simply to set 
to zero the lowest singular values. The condition number 
becomes: 




1 p_9A E 1 I 1 < I 1 1 — I 3 

10 100 1000 

rank 

Fig. 9. Typical set of singular values for a 36 segment ger- 
manium coaxial crystal (boxes: 5 mm cubic grid, dots: 2 mm 
cubic grid). The signal are generated using the MGS code. 



As can be seen the transformation matrix of a coax- 
ial detector cannot be characterized by a numerical rank. 
Germanium signal decomposition is an ill-posed problem. 
Hence, the value of r will have to be tuned in order to find 
the right balance between the precision of the solution vec- 
tor Gsoi and the amount of fluctuations on its components. 

We now turn to the implementation of this condition- 
ing improvement method. 
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3.3 SVD truncation 



The handling of the matrices in SVD truncation is illus- 
trated in fig JlO] which shows graphically the matrix prod- 
uct of eq. (I7|). The first step of the method consists in 
truncating the system, keeping only the relevant r largest 
singular values. The red hatched part of the W matrix 
is now full of zeros. Hence, being multiplied by zeros, the 
lower part of the *V and the right part of the U matri- 
ces can also be discarded (hatched blocks). The resulting 
reduced matrices are noted , and W^. As can be 
seen, this operation does not modify the size of the M 
matrix. The goal here is to improve its conditioning. 




■ 












m 







M 



Fig. 10. Graphical sketch of SVD truncation. The smallest 
singular values are replaced by zeros, thus the red hatched 
blocks of matrix W are made of zeros only. Therefore, the 
bottom part of the W*V matrix is also made of zeros. Being 
multiplied by zeros, the bottom part of matrix *V and the 
right part of matrix U play no role. SVD truncation consists 
in removing the hatched areas from the matrices. The resulting 
reduced matrices are U^, and ^V^. 



Doing this, we have, of course, introduced a systemat- 
ical bias on the solution but, at the same time, we have 
reduced the uncertainty on the solution. The optimum 
number of discarded singular values can be defined as the 
value corresponding to the minimum average square error 
on the hit location. This value can only be found empir- 
ically, varying r and calculating the solutions for a set 
of known locations. In fact, in most pulse shape analy- 
sis problems of germanium detectors, a large proportion 
of the singular values can be set to zero. Indeed, it has 
been observed that increasing the sample duration a lot 
has little effect on the precision of the hit location. The 
autocorrelation between successive samples is very strong 
(71/70 ^0.99 where 7/^ is the k^^ order autocorrelation 
coefficient), thus grouping the samples does not affect the 
information. The number of relevant samples being small, 
the number of relevant singular values is small as well. 



Another positive effect of setting to zero singular val- 
ues is that the number of lines in eq. ( 10 ) is lowered (see 



fig. 8k, which reduces drastically the computer time neces- 



sary to solve it. Indeed, as shown in fig. [10) the last lines of 
*V can be discarded from the matrices. This corresponds 
to a projection of the signal samples onto an optimum 
subspace spanned by the column- vectors of * that per- 
mits to keep a maximum information using a minimum 
number of components [3T] . 



3.4 Example of application 

In order to illustrate the method, we have applied it to 
determine the position of single hits in a given segment of 
an AGATA crystal. The basis signals as well as the test 
signals are given by the MGS simulation. The test signals 
are altered by adding different amounts of noise and by a 
3 ns standard deviation time shift jitter. A typical meta- 
signal is shown in fig. fl] The signal basis (M matrix) is 
made of 950 meta-signals of 468 samples corresponding to 
a 2 mm cubic grid covering the volume of the segment. 
Its condition number being of the order of 10^^, the direct 
solving of the system induces large uncertainties on the 
location of the hits. The distribution of the singular values 
shown in Fig 11 has no abrupt gap, thus the optimum 
number of discarded singular values has to be determined 
empirically. 
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Fig. 11. Distribution of the singular values for a single AGATA 
detector segment for a 2 mm cubic grid. 



The truncated system, eq. (10) and fig. 10, reads: 



^V'e = s' with s' = W'"^*U's. 



(13) 



This equation can be solved in several ways. One is the 
use of a fast iterative inversion algorithm. We have used 
NNLS. Another method consists in comparing the reduced 
test signal to all the signals of the truncated basis of re- 
duced signals (grid search on reduced signals). Here, 
the criterion for the best match is the residue [29 (the 
scalar product criterion [19 results in a faster, but slightly 
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less precise, algorithm). The results obtained with both 
methods and for different signal-to-noise ratios are shown 



in fig. 12 
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Fig. 12. Average error on the hit location as a function of 
the number of retained singular values, in the case of single 
hits at random positions in a given segment of an AGATA 
detector. Top panel: grid search; bottom panel: NNLS min- 
imization. The line styles correspond to different amount of 
noise. Bold: SNR = 42.4 dB, thin: SNR = 32.4 dB, dashed: 
SNR = 22.4 dB, dot-dashed: SNR = 12.4 dB, dotted: 
SNR = 2.4 dB. A 3 ns time jitter is also applied to the test 
signals. The vertical bar corresponds to 16 singular values. 



In every case, the errors are large when only a few 
singular values are kept since the remaining information 
is too small to allow a precise localization, and when the 
number of singular values is large since the condition num- 
ber is large. The optimum number of singular values de- 
pends slightly on the algorithm and on the signal-to-noise 
ratio. The resolution when 16 singular values are retained 
is close to the minimum, whatever the inversion method, 
the SNR, and for the time shift values that can be ex- 
pected from electronics. It is remarkable that the optimum 
number of singular values is only twice the number of sig- 
nals included in the meta-signals. For both algorithms, the 
computing time is proportional to the number of singular 
values. The SVD method, however, entails the calculation 



of the reduced meta-signal event by event. Nonetheless, 
keeping 16 singular values out of 468 reduces the compu- 
tation time by a factor of about 20. 

This analysis also shows that grid search is well adapted 
to the determination of single hits, since it is the fastest al- 
gorithn]^ The inversion method developed in the next sec- 
tion gives better results when the signals, resulting from 
several hits, overlap. 



4 Signal decomposition algorithm 

4.1 More than one hit in a single segment 

Depending on the type of application, the situations when 
a gamma interacts more than once in the same segment (or 
two gamma interact simultaneously in the same segment), 
are treated in different ways. For source location, or more 
generally, when the precision on the measurement is more 
important than the amount of analyzed photons, these 
events are simply discarded. However, even in this case, 
in order to be rejected, multi-hits have to be discriminated 
from single hits. 

In other cases, typically in nuclear physics experiments 
yielding high gamma-ray multiplicities, discarding events 
would introduce a crippling bias. Hence, whenever possi- 
ble, every hit has to be located. Grid search algorithms 
are not well adapted to the solving of multi-hits. Indeed, 
for single hits the computing time is proportional to the 
number of grid points, for double hits this time is propor- 
tional to the square of the number of grid points times the 
number of possible energy sharing between the two hits. 
Thus the computing time increases more rapidly than the 
exponent of the number of hits. This drawback is also true 
for most artificial intelligence techniques. Using the ma- 
trix formalism, as will be seen now, a faster algorithm can 
be developed. 



4.2 Locations of the multi-hits 



Signal decomposition consists in solving eq. (13) in or- 



der to determine the location and the deposited energy 
for each gamma interaction. In the multi-hit case, each 
hit appears in the solution vector as a cloud of adjacent 
(in the position space) non-zero components (fig. 13). If 
the clouds have strong overlaps then it is not possible to 
distinguish the different hits. This happens when the hits 
are too close or when they occur in a part of the segment 
where the sensitivity is low or when the energy deposits 
are low (see low SNR part of fig. [4|. When the clouds are 
separated, each hit location is estimated using eq. ([5| in 
which the sum runs only on the components of the cor- 
responding cloud. Both the cloud discrimination and the 
hit location estimation can be computed rapidly using the 



^ Only when the noise is low (SNR > 18 dB), the Grid 
Search solving of eq. Q may give a slightly better resolution 



than eq. (13) truncated to 16 singular values, but the comput- 



ing time is also 20 times longer. 
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mobile center algorithm. This method is illustrated in fig. 
[T3I The estimated locations of the hits are first sorted 
at random. Then each positive component of the solution 
matrix is associated to the closest of these centers and 
each center is replaced by the barycenter of its associated 
components. This procedure is repeated iteratively until 
the center locations are stable. Their final positions are 
the estimated locations of the hits. 



a) 



o 



b) 



4.3 Protocol for the whole crystal 

We now consider the whole array of segments. Instead of 
being built with the signals of the hit segment and of its 
neighbors, the meta-signals are obtained by the concate- 
nation of all the segment signals and the vector e gather 
the energy deposits on all the crystal grid voxels. This cor- 
responds to a large increase of the transformation matrix 
size. However, using SVD truncation, the resulting size is 
not excessive for on-line applications. 



III 




1 








hi 




1 








t 


: 





M 



Fig. 14. Graphical sketch of singular value decomposition af- 



Fig. 13. Illustration of the mobile center algorithm. The sur- 
faces of the dots are proportional to the corresponding com- 
ponent in the e vector; their colors correspond to the closest 
mobile center. The two mobile centers are represented with di- 
amonds, a) Initial (random) positions of the mobile centers, b) 
positions after one iteration, c) positions after two iterations 
(final positions) and estimated hit locations. 



Most of the time, the main difficulty is to determine 
the number of simultaneous hits in a segment. Many a 
priori algorithms (in the sense that the number of hits is 
determined first) have been tested but the performances 
are still not very satisfactory. An interesting strategy is 
proposed in reference [34,- The matrix formalism allows 
an a posteriori determination. First, the linear system is 
solved. The total cloud of non-zero components can then 
be analyzed. Several strategies are possible for searching 
for independent sub-clouds corresponding to the different 
hits. We have found that the analysis of the cloud inertia 
tensor (multi-hits correspond to larger moments) was the 
most robust method in the case of coaxial detectors. 

An important advantage of the singular value decom- 
position method is that it reduces the dispersion of the 
clouds. This property facilitates the separation of the clouds, 
thus the identification of the hits. Yet, whatever the cho- 
sen method, it seems difficult to discriminate more than 
two hits in a segment and the errors on the location and 
the energy sharing are small only if the locations are not 
too close and the energies are not too low. 



ter truncation (the red hatched blocks of fig. 10 are discarded). 
In typical events, not all the segments are hit simultaneously, 
thus, only the columns of M (dotted line) and * correspond- 
ing to the hit segments have to be kept for the location of the 
hits. Moreover, when given segments are hit, signals are seen 
only for these segments and their neighbors. Thus, only the 
corresponding lines of the M and matrices have to be used 
in the signal decomposition (green blocks). 



In a typical gamma event, only a small number of seg- 
ments are hit, thus it is not necessary to solve the linear 
system using the whole M and e matrices. Only the co- 
ordinates of the unknown vector corresponding to grid 
points inside the hits segments, and the corresponding 
columns of the basis signal matrix, have to be retained. 
In fig. 14, we consider a case where the hit segments cor- 
respond to the two columns indicated by dotted lines in 
matrix M. The interactions induce signals only in the hit 
segment and its neighbors. The useful part of the meta- 
signals are indicated by the green blocks in the M matrix. 
Thus, we introduce the *V^^* and U^^* matrices which are 
built only with the useful blocks. Similarly, the e^^* vector 
is composed only with the grid points belonging to the hit 
segments (green columns of * V^) and s^^* is composed by 
the signals of the hit segments and their neighbors (green 
lines of U^). 

As we have seen, the system to be solved is given by 
eq. (13) and fig.js] The R^^* matrix is calculated first (fig. 
15). The final linear system is shown in fig. [16] As can be 
seen, the size of the system has been drastically reduced. 

In order to use this method in on-line applications, it is 
important to pre-calculate as many matrices as possible. 
The singular value decomposition can be done off-line and 
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Fig. 15. Graphical sketch of the construction of the s^^^* re- 



duced signal corresponding to the right-hand side of eq. ( 13 ) 



I 

□ I 

t^-hit 



Fig. 16. Graphical sketch of eq. ([13|. The ^V^'* matrix is the 
concatenation of the useful (green) blocks of matrix ^V^. Only 
the components of vectors e corresponding to the grid points 
which are inside one of the hit segments, are retained. 



For high-rate on-line applications, the SVD of the whole 
transformation matrix is computed off-line. Each event 
has to be solved selecting only the relevant blocks of the 
matrices. 

In conclusion, whatever the algorithm used for signal 
decomposition, the signal basis should be analyzed and 
reduced using SVD in order to speed up the on-line com- 
putations and to improve the conditioning of the response 
function. For complex events, that is when the detected 
signal results from multiple gamma interactions, the ma- 
trix method gives very good results in terms of energy and 
position precision as well as of computing time. 



A Appendix A: Box diagram for matrix 
products 



the resulting and R matrices are memorized. Event 
by event, the ^ V^^* and R^^* matrices have to be extracted 
from the previous matrices and the reduced signal is ob- 
tained by multiplying the detected signal by the R^^* ma- 
trix. Finally, only the small * V^^* e^^* = system is to 
be solved. 

One advantage of this protocol is that the singular 
value decomposition has to be done only once, off-line. 
Thus, on an event-by-event basis, one has only to select 
the useful parts of the matrices and solve a small system 
inversion. 

A more rigorous, but more complex, procedure is first 
to select the useful blocks then to the calculate the sin- 
gular value decomposition. This protocol is developed in 
Appendix B. 



5 Conclusion 

The response function of High Purity Ge detectors has 
several characteristics that make signal decomposition dif- 
ficult. The relation between the pulse shapes and the lo- 
cations of the hits is not always bijective, the response 
function amplifies the signal noise to the hit location esti- 
mation, signal decomposition is an ill-posed problem and 
the size of the linear system to be solved is very large. 

These problems have been addressed using the matrix 
formalism. The first advantage of this method is to allow a 
mathematical analysis of the response function of the indi- 
vidual segments or of the whole germanium crystal. Using 
the SVD analysis, we have evaluated how signal uncer- 
tainties alter the precision on the estimated hit location. 
This decomposition also indicates the maximum number 
of grid points that a segment or the crystal can accommo- 
date. We then discuss the SVD truncation leading to the 
reduction of the size of the system, which permits both 
to reduce the computing time and to decrease the uncer- 
tainty on the hit location estimation. 



For a better readability, the matrix handling of this paper 
are illustrated by box diagrams. This Appendix shows how 
such diagrams are constructed. A matrix is represented 
by a box which height is proportional to the number of 
lines and which width is proportional to the number of 
columns. The product of two matrices is represented in 
the following way: 



B 


bij 

bnj 






Q>il . . . din 













A C 
Fig. 17. Graphical sketch representing the product A B = C. 



The Cij component is the sum-product of the line of 
A and the column of B pointing towards it. The vectors, 
assimilated to single column matrices, are simply repre- 
sented by a vertical line as in fig. [3] Some properties of 
matrix products can be easily visualized using box dia- 
grams. For example, as shown in fig. [Tsj if the last line of 
A is made of zeros, then the last line of C is also made 
of zeros. If the last line of B is made of zeros then the 
removals of the last column of A and of the last line of B 
do not modify C. 

The product of three matrices can be represented by 
two different diagrams (fig. 19). 
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Fig. 18. Illustration of matrix products properties using a box 
diagram. The red hatched blocks of A and B can be discarded. 
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Fig. 19. Diagrams for the product of three matrices: ABC 
D. 
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B [-'^ □ I 



S2 



MS 



Fig. 20. Splitting of the linear system into two independent 
smaller systems, a) Selection of the useful blocks of the initial 
M matrix (no block overlap as in fig. 
independent linear systems. 



14). b) Resulting two 



B Appendix B: 
segments 



Protocol for subsets of 



B.l Splitting of the linear system 

An event is often composed of groups of hits which happen 
in well separated zones of the detector, in that sense that 
the signals induced by each group do not overlap with the 
signals induced by the others. In that case, every group of 
hits can be treated independently, fig. |20] shows how the 
matrix system can be split into two separated systems. 



B.2 Inverse protocol 

A more rigorous procedure to solve the decomposition 
problem is first to select the useful blocks of the transform 
matrix then to calculate the singular value decomposition 
(fig. 21 ) for the selected blocks. In this case, the singular 
value decomposition has to be performed event-by-event, 
thus this methode is appropriate for off-line applications 
or when the acquisition rate is low. 

However, for a reduced number of segment combina- 
tions (the most probable ones) , it is possible to pre-calculate 
off-line the *V^^* and R^^* matrices resulting from the se- 
lection of the blocks followed by their singular value de- 
composition. This set of matrix pairs is saved in the com- 



puter memory so that it can be used on-line when the cor- 
responding segment combination is encountered. For the 
combinations that were not memorized, the previous pro- 
tocol is used. Table I shows that, due to computer memory 
limitation, this clustering method can be used for coaxial 
detectors only when each combination involves less than 
three or four adjacent hit segments. 



Table 1. Number of possible segment combinations as a func- 
tion of the number of hit segments in the case of 36 segment 
coaxial detectors. A combination gather segments which are 
close enough for their signals to overlap. 



number of 


number of 


hit segments 


combinations 


1 


36 


2 


315 


3 


3 146 


4 


22 951 


5 


137 957 
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Fig. 21. Inverse protocol for sub sets of the detector, a) Se- 
lection of the useful blocks of the initial matrix M. The blocks 
labeled by a are composed of zeros only, b) Building of the 
M^^* matrix from the selected blocks and SVD truncation, c) 
Reduced system to be solved (R^^* = W^^'"^ *U^^'). 
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